# coding=UTF-8

from scipy.ndimage import filters
from numpy import *

def compute_harris_response(im,sigma=3):
    '''计算Harris角点检测器的响应函数
    for each pixel
    '''

    imx=zeros(im.shape)
    #执行高斯滤波,参数(0,1)表示先对原图像作高斯核(零阶)卷积，再将图像与高斯核的一阶导做卷积
    filters.gaussian_filter(im,(sigma,sigma),(0,1),imx) #x方向求导
    imy=zeros(im.shape)
    filters.gaussian_filter(im,(sigma,sigma),(1,0),imy) #y方向求导

    #计算Harris矩阵的各元素
    Wxx=filters.gaussian_filter(imx*imx,sigma)
    Wxy=filters.gaussian_filter(imx*imy,sigma)
    Wyy=filters.gaussian_filter(imy*imy,sigma)

    #行列式与迹
    Wdet=Wxx*Wyy-Wxy**2
    Wtr=Wxx+Wyy

    return Wdet-0.05*Wtr*Wtr


def get_harris_points(harrisim, min_dist=10, threshold=0.01):
    """从一个Harris响应图像返回角点
    min_dst是隔离图像边界与角点的最小像素个数
    """
    #找到在阈值之上的最大的候选角点
    corner_threshold=harrisim.max()*threshold
    harrisim_t=(harrisim>corner_threshold)*1    #乘1将逻辑True，False转为1和0

    #得到候选点的
    coords=array(harrisim_t.nonzero()).T
    #以及对应的值
    candidate_values=[harrisim[c[0],c[1]] for c in coords]

    #对候选点排序，得到的是candidate_values里面的值的排序索引。
    index=argsort(candidate_values)
    #创建一个允许区域，在此区域内的角点才有效
    allowed_locations=zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist]=1

    #选择角点
    filtered_coords=[]
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] ==1:
            filtered_coords.append(coords[i])
            #此步骤保证在当前区域中只取一个有效角点
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
                 (coords[i,1]-min_dist):(coords[i,1]+min_dist) ]=0
    return filtered_coords

def plot_harris_points(image,filtered_coords):
    """绘制在图像中打到的角点。"""
    figure()
    gray()
    imshow(image)
    plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords],'r*')
    axis('off')
    show()


from pylab import *
from PIL import Image

im=array(Image.open(r'misc_pic\chessboard.jpg').convert('L'))
harrisim = compute_harris_response(im)
filtered_coords=get_harris_points(harrisim,3)
plot_harris_points(im,filtered_coords)